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ABSTRACT 

We study the effect of varying tlie equation of state on tlie formation of stellar 
clusters in turbulent molecular clouds, using three-dimensional, smoothed parti- 
cle hydrodynamics simulations. Our results show that the equation of state helps 
determine how strongly self-gravitating gas fragments. The degree of fragmenta- 
tion decreases with increasing polytropic exponent 7 in the range 0.2 < 7 < 1.4, 
although the total amount of mass accreted onto collapsed fragments appears to 
remain roughly constant through that range. Low values of 7 are expected to lead 
to the formation of dense clusters of low-mass stars, while 7 > 1 probably results 
in the formation of isolated and massive stars. Fragmentation and collapse ceases 
entirely for 7 > 1.4 as expected from analytic arguments. The mass spectrum 
of overdense gas clumps is roughly log-normal for non-self-gravitating turbulent 
gas, but changes to a power-law under the action of gravity. The spectrum of 
collapsed cores, on the other hand, remains log- normal for 7 < 1, but flattens 
markedly for 7 > 1. The density PDFs approach log-normal, with widths that 
decrease with increasing 7. Primordial gas may have effective 7 > 1, in which 
case these results could help explain why models of the formation of the first 
stars tend to produce isolated, massive objects. 



Subject headings: ISM: clouds — ISM: hydrodynamics — stars: formation — 
equation of state — turbulence 
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1. INTRODUCTION 

Stars form alone and in groups, associations, and clusters (see Pudritz 2002 and Ward- 
Thompson 2002 for recent reviews). Yet the origin of this diversity in the process of star 
formation remains poorly understood. Recent work suggests that stars form by gravitational 
collapse of shock-compressed density fluctuations generated by the supersonic turbulence 
generally observed in molecular clouds (e.g. Elmegreen 1993; Padoan & Nordlund 1999; 
Klessen, Heitsch, & Mac Low 2000; Ostriker, Stone, & Gammie 2001; see Mac Low & 
Klessen 2003 for a recent review). If so, then star formation is controlled by turbulent 
fragmentation and dynamical collapse of individual, Jeans-unstable, protostellar clumps. 

Fragmentation of gas clouds has been studied for more than a century, yet the process 
remains poorly understood. There have been analytical models (Jeans 1902; Low & Lynden- 
Bell 1976), and numerical investigations of the effects of various physical processes on the 
collapse, such as geometry and rotation of the clouds (see Bonnell & Bastien 1993 for a 
review), and magnetic fields (Galli et al. 2001). Recently, the effects of turbulence have been 
studied extensively in a series of three-dimensional simulations using both grid-based, Eule- 
rian and particle-based, Lagrangian hydrodynamics (Klessen, Burkert, & Bate 1998; Klessen 
& Burkert 2000, hereafter Paper I; Klessen, Heitsch, & Mac Low 2000, Paper II; Heitsch, 
Mac Low, & Klessen 2001, Paper HI; and Klessen 2001, Paper IV). These papers show that 
turbulence strongly influences the fragmentation of molecular clouds. Fast, clustered collapse 
and star formation occurs in regions with turbulence insufficient to support against gravita- 
tional collapse, while slow, scattered star formation results from strong turbulent support. 
However, these results are based on isothermal models with polytropic exponent 7 = 1 only. 
The role of the equation of state (EOS), which is essential in understanding the physical 
structure and stability of the turbulent clouds, remained unexplored. 

The balance of heating and cooling in a molecular cloud can be approximately described 
by a polytropic EOS 

P = Kp\ (1) 

where X is a constant, and P, p and 7 are thermal pressure, gas density and polytropic 
exponent, respectively. The density structure generated by supersonic motions in highly 
compressible turbulence depends on the EOS. The stiffness of the EOS can largely con- 
trol the density probability distribution function (PDF) in strongly compressible turbulence 
(Scalo 1998; Passot & Vazquez- Semadeni 1998; Vazquez-Semadeni & Garcia 2001). In par- 
ticular, one-dimensional simulations by Passot & Vazquez-Semadeni (1998) show the shape 
of the density PDF varying with 7. They found that the PDF is log-normal when 7 = 1 
(isothermal), and develops a power-law tail for 7 7^ 1. Spaans & Silk (2000) showed that 
radiatively cooling gas can be described by a piecewise polytropic EOS, in which the poly- 
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tropic exponent 7 changes with gas density p. They conducted a detailed analysis including 
chemistry, thermal balance and radiative transfer, and found that in the interstellar medium 
the effective polytropic exponent 7 has a range of 0.2 < 7 < 1.4. They further predicted that 
the initial mass function (IMF) of protostellar cores would vary with 7, yielding a peaked 
IMF for 7 > 1 and a power-law function for 7 < 1. 

In this paper, we perform three-dimensional simulations of sclf-gravitating, supersonic 
turbulence over the range of 7 identified by Spaans & Silk (2000) to investigate the effect of 
varying the EOS on the dynamical evolution and fragmentation behavior of turbulent clouds 
and the subsequent formation of protostellar clusters. To quantify the dependency on 7, we 
consider a true polytropic EOS, with 7 being strictly constant in each simulation regardless 
of the density, rather than the piecewise polytrope used by Spaans & Silk (2000) and earlier 
workers. 

In §2 we describe our computational methods; in §3 we give results on fragmentation, 
the resulting mass distribution of protostellar cores, and the gas density PDF; and in §4 we 
conclude and speculate about the implications or our work for the formation of star clusters 
and primordial star formation. 



In analytical work on the stability of a self-gravitating, isothermal medium, Jeans (1902) 
derived a relation between the oscillation frequency uj and the wavenumber k of small per- 
turbations. 



where is the isothermal sound speed, G the gravitational constant, and p is the gas density. 
The medium is unstable to fragmentation at all wavelengths greater than a critical length 
Aj = 27r//cj, or equivalently all masses exceeding the Jeans mass. 



2. 



METHODS 



u'^ - c'^e + AttGp = 



(2) 




(3) 



will collapse under their own weight. Note that we use a cubic definition of Mj. 



In a polytropic cloud with EOS given by equation 1, the sound speed is 




(4) 



so the Jeans mass will be given as a function of 7 and p, 




(5) 
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The specific internal energy is it = Kp^ ^ /{lad ~ 1); where ^ad is the adiabatic index, with 
lad = 5/3 for ideal gas. For isothermal gas, 7 = 1, and K — c^^. 

We carry out simulations to determine when and how fragmentation occurs. We use a 
smoothed particle hydrodynamics (SPH) code (Benz 1990) in order to resolve several orders 
of magnitude in density during the collapse of shock-compressed regions. Our version of 
the code includes periodic boundary conditions (Klessen 1997) and can replace high-density 
cores with sink particles (Bate, Bonnell, & Price 1995). Sink particles accrete surrounding 
gas particles while conserving mass and momentum, but they only interact gravitationally. 
They prevent time step from becoming prohibitively short in very dense regions. This allows 
us to follow the dynamical evolution of the system over many free-fall times. 

We set up a periodic cube with side length L = 2 and total mass of unity such that 
initial density po = 0.125, and sound speed Cg = 0.1. The treatment of the polytropic 
EOS follows Bonnell & Bastien (1991), where we set K = {Rg/ fj,)TopQ~'^ . We choose values 
consistent with the sound speed of Rg — 2/3, A* = 1, and Tq = 0.015 for the gas constant, 
mean molecular weight, and initial temperature of the cloud. Uniform turbulence is driven 
with the method described by Mac Low et al. (1998) and Mac Low (1999), adding energy 
over a small range of wavenumbers k. Driving in two wavenumber ranges is considered here, 
1 < k < 2 and 7 < k < 8, corresponding to models Blh and B3 in Paper 11. In each case, 
the driving strength is chosen to ensure the same rms velocity despite the different driving 
wavenumbers, so that the average turbulent Jeans mass {Mj)turb = 3.2 and Jeans Number 
Nj = 64, are the same. In a self-gravitating medium, the maximum resolvable density is 
determined by the requirement that the local thermal Jeans mass be resolved by the SPH 
kernel (Bate & Burkert 1997); our models all have the same resolution limit. 

The models presented here are computed in normalized units. If scaled to mean densities 
of n(H2) ~ lO^cm^^, and a temperature of 11. 4K (i.e. an isothermal sound speed Cg = 
0.2kms~^), values appropriate for a dark cloud like Taurus (Hartmann 1998), then our 
simulation cube holds a mass of M 4 x 10^ M© and has a size of L 9 pc. For details on 
the scahng relations, see Paper II. 

We computed models with 0.2 < 7 < 1.4, covering the range suggested for interstellar 
gas by Spaans & Silk (2000). We vary 7 in steps of 0.1 in otherwise identical simulations 
for both driving wavenumbers. Each simulation starts with a uniform density po — 0.125. 
Driving begins immediately, while self-gravity is turned on at t=2.0, after the turbulence is 
fully estabhshed (see Paper II). Our models use 200 000 particles. Twenty-six simulations 
were performed simultaneously using a serial version of the code on single 1 GHz Pentium HI 
processors of the Parallel Computing Facility of the American Museum of Natural History 
(AMNH). Each computation took about 6 months. 
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3. RESULTS 

3.1. Turbulent Fragmentation 

Figure 1 and Figure 2 show the density distribution of the gas. To make the figures, a 
256^ grid was filled with densities properly computed from the SPH kernel, and the result is 
displayed in three-dimensional projection {top row) and in a slice through maximum density 
one free-fall time after self-gravity is turned on for selected 7 and both driving wavenumbers. 
As found in previous studies (e.g. Mac Low 1999), the driving wavenumbcr strongly influences 
the density distribution. Driving with k = 1-2 produces strong filamentary structure, while 
the density distribution for k = 7-8 remains more uniform at large scale. 

We find that the value of 7 is as important as the driving in determining collapse 
behavior. As 7 increases, the number of collapsed cores replaced by sink particles decreases, 
and the cores cluster less. For the same 7, driving with k = 1-2 produces both far more, 
and more clustered dense cores than driving with k = 7-8. The dense cores tend to collapse 
predominantly in filaments or at intersections of filaments. At 7 > 1.1, no fragmentation 
occurs by the time shown for driving with k = 7-8. 

The isothermal case, 7 = 1.0, agrees well with previous results in Papers 11 and 111. 
Our current investigation agrees with the results in Paper I-III, that collapse tends to form 
clusters with high efficiency in regions with weak turbulence, while in regions with strong 
turbulence, sparse, slow collapse occurs. We now add the additional criterion that collapse 
and fragmentation depends strongly on 7. 

In regions with 7 < 1, fragmentation occurs earlier and more frequently, while in regions 
where 7 > 1, fragmentation is retarded and less frequent. 

Figure 3 compares the number of collapsed cores at different 7 for models with different 
driving wavenumbers. The rate at which new protostellar cores form differs for different 7. 
Models with low 7 form cores quickly, while models with large 7 form new cores more rarely. 
Again we see more cores for driving with k — 1-2 than for k = 7-8. 

Rouleau & Bastien (1990) found an empirical scaling relation oc 7^ of the critical 
Jeans number Jc = M/Mj for collapse for elongated clouds. The Jeans number should be 
related to the number of fragments. However, with higher resolution and more statistics for 
a wide range of 7 in our simulations, we find something more like the inverse, but strongly 
dependent on the details of the driving. Fig 4 shows the relation between the value of 7 
and the number of fragments at one free-fall time. The number of fragments drops quickly 
as 7 increases for driving with k = 1-2, and more slowly for driving with k — 7-8, which 
produces fewer fragments to begin with. 
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Figure 5 shows the accretion histories (the time evolution of the combined mass fraction 
of all protostellar cores) for all the individual 7 cases and for both driving wavelengths. Note 
that much less mass is accreted when driving with k = 7-8. We focus on the k = 1-2 model. 
The lower 7 is, the earlier fragmentation occurs and the larger the number of fragments that 
form. It also takes longer for the clouds with high 7 to fragment. The slope of the accretion 
curves are roughly the same, but the number of sinks is very different for different 7 at the 
same time, which suggests that the average core mass is different for different 7. A low-7 
environment produces a large number of collapsed cores of relatively low mass, while 7 > 1 
results in the formation of considerably fewer, but more massive, cores. Since the overall 
mass growth rate of collapsing cores is closely related to the expected star formation rate, our 
result suggests that the overall star formation rate is higher in low-7 clouds than in high-7 
clouds. It remains unclear what the final amount of mass in stars will be, as higher mass 
stars may drive surrounding gas away more quickly (see also Vazquez-Semadeni, Ballesteros- 
Paredes, & Klessen 2003). 

We do not find any signs of collapse in the model with 7 = 1.4 even after several free-fall 
times. This can be explained with a Jeans mass analysis giving the critical condition for 
gravitational collapse. From Equation [5], we have 

This imphes that for 7 < 4/3 the Jeans mass Mj decreases as density increases during 
collapse; for 7 = 4/3, Mj remains constant; while if 7 > 4/3, then Mj increases with 
density. Thus, any collapse will choke itself off when 7 > 4/3, as found in our numerical 
results. 

Fragmentation is a complex process depending on the local conditions in the cloud, 

so the above analysis for an isolated spherical perturbation is an approximation. Indeed, 
fragmentation already ceases for 7 > 1.1 in the model driven with A; = 7 — 8, suggesting the 
importance of the details of the turbulence in determining collapse behavior. 



3.2. Mass Spectrum 

Observations suggest that the mass distribution of gas clumps in molecular clouds follows 
a power law, dN/dM = M'^, with typical values of the exponent being z/ ^ —1.5 (e.g. 
Williams, Blitz, & McKee 2000). Salpeter (1955) derived a power law for the high mass 
stellar IMF in the same notation with 1/ = —2.35 (see also Scalo 1998; Kroupa 2002). 

We measure clump-mass spectra in our models using a clump-finding algorithm similar 
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to the one described by Williams, De Geus, & Blitz (1994) but working on all three spatial 
coordinates and adapted to make use of the SPH kernel smoothing procedure (for details 
see Appendix A in Paper I). Figure 6 shows the spectra of gas clumps and collapsed cores 
for models with different values of 7, driven with k = 1-2. Three different evolutionary 
phases are shown: after turbulence has been fully established, but before self-gravity has 
been turned on, at t = 2; and when the fraction of mass accumulated in collapsed cores (sink 
particles) has reached ^ 20% and 40%. 

The mass distribution of both clumps and collapsed cores changes with 7, with the 
effect being most pronounced for the cores. In low-7 models, the core mass spectrum at the 
high-mass end is roughly log-normal. As 7 increases, fewer but more massive cores emerge. 
When 7 > 1.0, the distribution is dominated by a few high mass cores, and the spectrum 
tends to flatten out. It is no longer fit by either a log- normal or a power-law. The clump 
mass spectra, on the other hand, do show power-law behavior on the high mass side, even 
for 7 > 1.0. 

Our results suggest that massive stars can form in small groups or even in isolation in 
gas with 7 > 1.0. Spaans & Silk (2000) suggested that a stiff EOS with 7 > 1.0 should 
lead to a peaked IMF, biased toward massive stars, while an EOS with 7 < 1.0 results in a 
power-law IMF, in general agreement with our simulations. 



3.3. Density Probability Distribution Function 

Figure 7 shows mass and volume weighted PDFs of gas density and before self- 
gravity has been turned on for selected values of 7 and both driving wavenumbers. We also 
show Gaussian fits to the PDFs together with their error e — '^\p — Pfit\/P: summed over 
all points above 10% of the peak values. 

The mass- weighted density PDF Pm{p) is calculated directly from the local density 
associated with each SPH particle, as described by Klessen (2000). The volume- weighted 
density PDF Pv{p) is calculated on a cubic grid derived from the SPH density field by taking 
a cube of 256^ cells and using the SPH kernel smoothing procedure to find the density 
at the center of each cell. Both Pm{p) and Pv{p) are normalized so that JT^Pm{p)dp — 
n^Pv{p)dp=l. 

The density behind shocks in supersonic turbulence depends on the compressibility of 
the gas, 

dp(x-P^-^dP, (7) 
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So turbulent density fluctuations should increase as 7 decreases. Indeed, we see that in 
Figure 7, the width of Pm and p„ increases as 7 decreases. 

In models driven with wavenumber A; = 7 — 8, there are super-Gaussian tails in the 
low-density end in the volume-weighted PDFs for all 7 and in the mass-weighted PDF for 
7 = 0.2. One-dimensional simulations by Passot Sz Vazquez- Semadeni (1998) showed that 
the volume- weighted density PDF of supersonic turbulent gas displays a power-law tail at 
high densities for < 7 < 1, becomes log- normal for 7 = 1 (isothermal), and develops a 
power-law tail at low density for 7 > 1, which is only in partial agreement with our results. 
The discrepancy may be due to the low resolution in low-density regions in our simulations, 
the difference between one- and three-dimensional simulations, or different driving. 

A full investigation of the density PDFs under various conditions, with different codes, 
different dimensions (three-dimensional vs. one-dimensional), different resolution, different 
turbulence driving, different initial Mach number, and with and without magnetic field is 
beyond the scope of this paper; these issues will be addressed elsewhere. 



3.3.1. Comparison Between Mass and Volume Weighting 

In order to compare the mass and volume weighted PDFs, we will find it useful to 
define s = ln(p/po); the mass and volume weighted PDFs of s are then Pm{s) and Pv{s), 
respectively. We can relate p{s) to p{p) if we remember that for any montonic function y{x), 
it can be shown that the PDF \p{y)dy\ = \p{x)dx\. Therefore, 

Pm(p) ^Pm(s)^ = , (8) 

dp p 

and 

P.(P)=P.(«)^ = ^- (9) 
dp p 

Since Pm(p) oc dM/dp, andp^(p) oc dV/dp (Ostriker, Stone, & Gammie 2001), it follows 

that 

dMdV 

Pm{P) oc -^-^ ppvip). (10) 

Then we can relate Pm{s) to Pv{s) 

Pm{s) = PPmip) = Cp^p^{p) = Ce'p^{s) . (11) 



where C is a normalization constant. 



-9- 



If Pv{s) follows a normal distribution, 

where Sy is the average (in volume) value of s, and a is the dispersion, then 

= ^exp f-^i^%±^l exp f-Ii^%t^ I . (13) 



/2^a V 2(72 y " \ 2(72 

Since J^^Pm{s)ds — 1, the normalization 



^^ = expr^?-4^±^) . (14) 



2(72 



Equation 13 can thus be rewritten as: 



pU«) = ^exp(-^^^) • 

where Sj^^ = s^, + cr^. This is also a normal distribution, with the same dispersion a as that 
of pv but shifted average value Sm- Using the volume PDF normalization J p • py{s)ds — 
J e^py{s)ds = 1, it can be seen 

as noted, for example, by Ostriker, Stone, & Gammie (2001). This holds only if the density 
PDF is log-normal. (Note that in this case C — 1.) 

This derivation gives several properties of the density PDFs p„i and p^j. If one of them is 
log-normal, so is the other; and if they arc log-normal, the width a of both profiles should be 
equal (Eqs. [12] and [15]), and the peaks of the PDFs Pm{s) andp„(s) should lie symmetrically 
around zero (Eq. [16]). 

In Figure 7 the volume-weighted PDFs lie to the left of the mass-weighted PDFs, with 
—Sy, as predicted. The decent Gaussian fits to the peaks of the PDFs show that the 
assumption of log-normal behavior is not bad. Three-dimensional grid-based simulations of 
decaying turbulence by Ostriker, Stone, & Gammie (2001) show similar Gaussian PDFs. 

A detailed examination of the moments of the density PDFs, shown in Fig 8, shows that 
these PDFs are, in fact, not perfect log-normals. The first moment (mean) shows that p^ 
and py of models with both driving wavenumbers A; = 1 — 2 and k — 7 — 8 are not exactly 
symmetrically distributed. This could be due to turbulent intermittency in our small boxes 
producing deviations from a purely Gaussian distribution. The Gaussian behavior expected 
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from a well sampled distribution is more visible in models driven with A; = 7 — 8 compared 
to A; = 1 — 2. In fact, Pm and py of k = 7 — 8 are more equal in width, as seen from the 
second moment (variance). This is a sign of better sampling - more modes contribute to 
the overall velocity field, thus the central limit theorem is more appropriately applied to the 
A; = 7 - 8 than to the A; = 1 - 2 case (Klessen 2000). 

3. 3. 2. Effect of Self- Gravity 

Figure 9 shows the evolution of the mass-weighted gas density PDF pm of the model 
driven with A; = 1 — 2 for different values of 7. The same three evolutionary phases as 
defined in §3.2 are shown: the initial distribution, and times at which dense cores have 
accreted M^, ^ 20% and M^, ^ 40% of the mass. These PDFs do not include gas accreted 
into the cores. We again show the best-fit Gaussian for the part of each PDF above 10% of 
peak. 

During the dynamical evolution, the PDFs in the non-isothermal cases develop pro- 
nounced positive deviations from the the Gaussian fit, while the PDF in the isothermal 
case remains well fit by a Gaussian (of increasing width), consistent with the results in 
Klessen (2000). The high-density tails gradually diminish as 7 increases from 0.2 to 1.0, 
vanishing almost entirely at 7 = 1.0. At 7 > 1.0, though, the tails develop again and get 
stronger as 7 continues to increase. We do not fully understand why the tails only appear 
in non-isothermal cases. 

The first four moments of the density PDFs for each of the three evolutionary phases 
discussed previously are shown in Figure 10. We can see that the moments vary as collapse 
proceeds, demonstrating that self-gravity plays an important role in shaping the density 
PDF. 

The plot of the skewness shows that the PDFs of turbulence without self-gravity are 
always skewed to low densities, regardless of 7. Shock compression comes at a cost: the gas 
swept up in shock fronts has to come from somewhere, and those regions now contain only 
low-density gas. As the filling factor of shocks is low (Smith, Mac Low, & Heitsch 2000) 
most of the volume is occupied by these low-density regions resulting in PDFs with negative 
skewness. The situation changes, however, when self-gravity comes into play. Gravitational 
contraction produces high densities in addition to the turbulent shock compression, the 
PDF's in the later stages of evolution are thus strongly biased towards high-density gas and 
the skewness becomes positive. 
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4. SUMMARY AND DISCUSSION 

We have investigated the effect of varying the equation of state (EOS) on the formation 
of stellar clusters in turbulent molecular clouds using a set of three-dimensional numerical 
simulations. Our results show that the EOS plays an important role in the fragmentation 
of the clouds, the determination of the initial mass function (IMF) of the protostellar cores, 
and the shape of the gas density PDFs. 

The ability of interstellar gas clouds to fragment under the action of self-gravity decreases 
with increasing polytropic exponent 7 in the range 0.2 < 7 < 1.4 relevant for Galactic 
molecular clouds. The total amount of mass in collapsed cores, however, appears to remain 
roughly constant through that range. At least half of this material is expected to accrete 
onto the protostellar system at the center of each collapsing core (e.g. Klessen 2001). Small 
values of 7 thiis lead to the formation of dense clusters of predominantly low-mass stars, 
while 7 > 1 results in the formation of isolated and massive stars. Fragmentation and 
collapse ceases entirely for 7 > 1.4 as expected from analytic arguments. 

The mass spectra of both clumps and collapsed cores change with 7, with that effect 
being most pronounced for the cores. For 7 < 1.0, their mass spectrum appears roughly 
log-normal. As 7 increases, fewer but more massive objects emerge. When 7 > 1.0, the 
distribution is dominated by high mass objects only, and the spectrum tends to flatten out. 
It is fit by neither a log-normal nor a power-law. The clump mass spectra, on the other 
hand, do show power-law behavior at the high mass end for all 7. 

The density PDF changes with 7 as well. The width of the profile decreases as 7 
increases, as expected from analytic arguments about shock compression. Both mass- and 
volume-weighted PDFs show imperfect log-normal distributions. Self-gravity helps shape 
the PDF. During dynamical collapse, the PDF develops a high-density tail, which is the 
imprint of local collapse of high density regions. 

In a polytropic gas cloud, any temperature fluctuation corresponds to a density fluctua- 
tion, A logiQ(T) ~ A logi^Q(p)(l — 7) —0.6. We note that from the gas density PDFs in Figure 
7, the maximum width of the profiles is about 12, which corresponds to Alogio(p) — 5. So 
for the extreme case 7 = 0.2, AlogiQ(T') ~ 3.5. For a cloud with an average temperature of 
about 10 K the maximum temperature of the low-density gas therefore could reach values 
up to 10^ K, if 7 is extremely low. This still falls in the temperature range observed in star 
forming regions. 

The principal motivation of this paper is to investigate how the fragmentation behavior 
of a turbulent cloud and the subsequent formation of stars are affected by physical state 
of the cloud, as described by the EOS. Since 7 reflects the balance between heating and 
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cooling, 7 is expected bo be low in regions of strong cooling, while 7 may be large where 
the cooling mechanism is insufficient (Spaans & Silk 2000). Unfortunately, the true chemical 
state and composition of interstellar gas is difficult to assess, so exact values of 7 are difficult 
to estimate observationally. 

The possibihty that isolated, massive stars form when 7 > 1 is of interest, not only 
because ionizing radiation, stellar winds, and supernovae from massive stars play an im- 
portant role in the interstellar medium, but also because massive stars are rarely observed 
in isolation. Rather, they are usually found as members of rich stellar clusters. Recently, 
however, Lamers et al. (2002) has reported observations of massive stars found in isolation 
or associated only with very small groups of lower- mass stars in the bulge of M51. Massey 
(2002) also reported finding apparently isolated, massive, field stars in both the Large and 
Small Magellanic Clouds. Our simulations suggest that molecular cloud regions forming 
isolated, high- mass stars have 7 > 1. In this case turbulent cloud fragmentation tends to 
produce isolated, massive collapsing cores that will feed onto isolated high-mass stars. Yorke 
& Sonnhalter (2002) demonstrated that high-mass stars can form in collapsing gas clumps 
via disk accretion very much the same way as low mass stars are believed to do. There 
is thus no need for a dense cluster environment where high-mass stars could build up by 
collisions of lower- mass stars (as suggested by Bonnell, Bate, & Zinnecker 1998). 

A theory of galactic-scale star formation partially based on the universality of log-normal 
density PDFs has been proposed by Elmegreen (2002). Our results do show a log-normal 
form at least for the peak of the PDFs, containing most of the gas. However, the dispersion 
of the PDF varies substantially as the EOS changes (Fig. 8), so that the fraction of gas 
available for star formation would depend significantly on the EOS in such a theory. This 
raises the question of how it could still be able to explain a global Schmidt Law across many 
different galaxies with different metallicities, radiation fields, interstellar pressures, and so 
forth. 

High resolution simulations by Abel, Bryan, & Norman (2000, 2002) of the formation of 
the first star suggest that fragmentation during the collapse of metal-free pregalactic halos 
is rather inefficient, resulting in the formation of single, massive stars rather than clusters 
of lower-mass objects. In the absence of metals, inefficient cooling may result in high 7, 
perhaps helping to explain these results. 

We do note, however, that simulations by Bromm, Coppi, & Larson (1999, 2002) show 
greater multiplicity. They use similar chemistry, but examine more massive, more isolated 
halos, with substantially larger limiting mass resolution. The EOS will have less effect on 
these larger scales, which are more dominated by gravity. If the metal-free gas indeed shows 
a high effective 7, our models would suggest that their collapsing regions will show little 
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further fragmentation if followed down to stellar mass scales. 
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log,o p 

Fig. 1. — Density distribution of the gas shown on a grid in three-dimensional projection 
{top row) and in a shce through maximum density, for driving with k — 1-2 at 1 t// after 
self-gravity is turned on. Values of 7 = 0.2 {left column), 7 = 1.0 {middle column), and 
7 = 1.3 {right column) are shown. The color bar for the density scale of the slices is at the 
bottom. Note the very high densities in these images apply to the collapsed cores only. The 
distribution of gas density is shown in Fig. 7. 
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Fig. 2. — Density distribution of tlie gas siiown on a grid in three-dimensional projection [top 
row) and in a slice through maximum density, for selected 7 and for driving with k — 7-8 
at 1 Tff after self-gravity is turned on. Values of 7 = 0.2 {left column), and 7 = 1.0 {right 
column) are shown. The color bar for the density scale of the slices is at the bottom. Note 
the very high densities in these images apply to the collapsed cores only. 
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Fig. 3. — Comparison of number of collapsed cores for models with different polytropic 
exponent 7, for driving with wavenumber k — 1-2 {left panel) and k — 7-8 {right panel). 
Note gravity was "turned on" at t = 2.0 ~ 1.33t//. The simulations of low-7 cases (0.2 and 
0.7 in this plot) of model k — 1-2 were terminated after a few t// due to prohibitively small 
time step. 
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Fig. 4. — Relation of the value of 7 to the number of collapsed cores one free-fall time after 
self-gravity is turned on. 
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Fig. 5. — Comparison of accretion rate of collapsing cores for turbulence driven with 
wavenumber k — 1-2 {left panel) and k — 7-8 {right panel). Note the different scales 
in the two plots. Gravity is turned on at t = 2.0 ~ 1.33r//. The simulations of low-7 cases 
(0.2 and 0.7 in this plot) of the model driven with k — 1-2 were terminated after a few Tff 
due to the prohibitively small time step. 
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Fig. 6. — Mass spectra of gas clumps {open histogram), and of collapsed cores {filled his- 
togram) for models with different 7 and driving with k — 1-2. The fraction of mass accreted 
into collapsed cores (sink particles) is given for the right two columns. The vertical line shows 
the SPH resolution limit. Also shown are two power-law spectra with v — —1.5 (dashed-line) 
and v = —2.33 (dotted line). 
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Fig. 7. — Mass-weighted {thick solid-line) and volume-weighted {thick dashed-line) gas den- 
sity PDFs Pm and of fully-developed turbulence prior to turning on self-gravity. Models 
with driving wavenumber k — 1 — 2 {left column) and A; = 7 — 8 {right column) are shown. 
Gaussian fits are shown with thin dotted and thin dot-dashed lines correspondingly, and the 
errors eu and ey summed over the parts of the curves above 10% of the peak values, are 
quoted. Note that s — ln(p/po)- 
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Fig. 8. — First moment (mean) /i (top) and second moment (variance) cr^ (bottom), of both 
Pm (filled) and py (open) as functions of 7, for models driven with wavenumbers k — 1 — 2 
(circle) and A; = 7 — 8 (square), where s — ln(p/po)- 
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Fig. 9. — Mass-weighted gas density PDFs Pm of models driven with A; = 1 — 2 for selected 
7. Only gas not accreted into sink particles is included. Three evolutionary phases are 
shown: the initial state {thick solid-line), 20% {thick dashed-line) , and 40% 

{thick dot-dashed line). The thin lines (solid, dashed, dot-dashed) are Gaussian fits to the 
corresponding PDFs above 10% of peak value. Again s — ln(p/po)- Note the change of the 
high-density tails with 7. 
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Fig. 10. — The first four moments of tliree evolutionary phases of the mass-weighted density 
PDF (corresponding to Figure 9)of the model driven with A; = 1 — 2, as functions of 7. , 
Shown are the mean /i (top left), variance (top right), the skewness S (bottom left), and 
the kurtosis of s = ln(p/ po). We define the fourth moment k, with a value of 3 subtracted, 
so that for a Gaussian k — 0, while an exponential has k — 3. 



